!
! Package for data conversion between Citizen and Julian. 
!
! 1996-12-23 F.Hroch
! 1996 to f90
! 2000 some change according to Numerical Recipes
!
! $Id$
!
! --------------------------------------------------------------------

Module trajd

  integer, parameter :: double = selected_real_kind(15)

contains

  !----------------------------------------------------------------------

  function datjd (year,month,day)

    ! Compute Julian date from input citizen year, month and day.
    ! Tested for all year

    real(double), parameter :: break = 15 + 31*(10 + 12*1582)

    real(double), intent( in ) :: year,month,day
    real(double) ::  datjd,y,m,d,a

    if( abs(year) < epsilon(year) ) stop "datjd: There is no Year Zero." 

    y = year
    if( y < 0.0 ) y = y + 1

    if( month > 2.0 )then
       m = month + 1.0
       d = day
    else
       y = y - 1
       m = month + 13.0
       d = day
    endif
    datjd = int(365.25*y) + int(30.6001*m) + d + 1720994.5
    if( d + 31*(m + 12*y) >= break )then
       a = int(y/100.0)
       datjd = datjd + 2.0 - a + int(a/4.0)
    endif
   
  end function datjd 

  !----------------------------------------------------------------------

  subroutine jdat(jd,year,month,day)

    ! compute citizen date: year, month and day from input julian date.
    ! only for jd>0! tested for all year except 1582-10-07..15.
    !
    ! WARNING! possibly work incorrectly for a negative year

    real(double), parameter :: break = 2299163.0

    real(double), intent( in ) ::  jd
    real(double), intent( out ) ::  year,month,day
    real(double) :: alpha,a,b,c,d,e
    
    if( jd >= break ) then
       alpha = int( ((jd - 1867216.0) - 0.25) / 36524.25)
       a = jd + 1.0 + alpha - int( alpha / 4.0 )
    else
       a = jd
    endif
    b = a + 1524.0
    c = int(6680.0+((b - 2439870.0) - 122.1)/365.25)
    d = 365.0*c + int(c/4.0)
    e = int((b - d)/30.6001)
    
    day = b - d - int(30.6001*e)
    month = e - 1.0
    if( month > 12.0 ) month = month - 12.0
    year = c - 4715.0
    if( month > 2.0 ) year = year - 1.0
    if( year <= 0.0 ) year = year - 1.0 
    if( year < 0.0 )then
       day = day - 0.5
    else
       day = day + 0.5
    endif
    
  end subroutine         

end Module trajd




